fish_traits_sum <- fish_traits %>%mutate(across(where(is.numeric), round, 2))## Display the table ----htmltools::tagList(DT::datatable(fish_traits_sum))
traits correlation
Code
M <-cor(numeric_traits[, c(-1)])ggcorrplot::ggcorrplot(M, hc.order =TRUE, type ="lower",lab =TRUE, tl.cex =9, lab_size =3)
Code
# list of species sp_names <-c(rownames(fish_traits))# taxonomic_familiestaxonomic_families <- sp_names %>%as.data.frame() %>%`colnames<-`("species") %>%mutate(family =case_when( species %in%c("Benthosema_glaciale","Ceratoscopelus_maderensis","Diaphus_metopoclampus","Lampanyctus_ater","Lampanyctus_crocodilus","Lampanyctus_macdonaldi","Lobianchia_gemellarii","Myctophum_punctatum","Notoscopelus_bolini","Notoscopelus_kroyeri","Bolinichthys_supralateralis" ) ~"Myctophidae", species %in%c("Borostomias_antarcticus","Chauliodus_sloani","Malacosteus_niger","Melanostomias_bartonbeani","Stomias_boa" ) ~"Stomiidae", species %in%c("Holtbyrnia_anomala","Holtbyrnia_macrops","Maulisia_argipalla","Maulisia_mauli","Maulisia_microlepis","Normichthys_operosus","Searsia_koefoedi","Sagamichthys_schnakenbecki" ) ~"Platytroctidae", species %in%c("Sigmops_bathyphilus","Gonostoma_elongatum") ~"Gonostomatidae", species %in%c("Argyropelecus_hemigymnus","Maurolicus_muelleri","Argyropelecus_olfersii" ) ~"Sternoptychidae", species =="Anoplogaster_cornuta"~"Anoplogastridae", species %in%c("Arctozenus_risso", "Paralepis_coregonoides") ~"Paralepididae", species =="Bathylagus_euryops"~"Bathylagidae", species =="Cyclothone"~"Gonostomatidae", species =="Derichthys_serpentinus"~"Derichthyidae", species =="Eurypharynx_pelecanoides"~"Eurypharyngidae", species =="Evermannella_balbo"~"Evermannellidae", species =="Lestidiops_sphyrenoides"~"Lestidiidae", species =="Melanostigma_atlanticum"~"Zoarcidae", species %in%c("Photostylus_pycnopterus","Xenodermichthys_copei") ~"Alepocephalidae", species =="Serrivomer_beanii"~"Serrivomeridae" ) )%>%mutate(order =case_when( family =="Myctophidae"~"Myctophiformes", family %in%c("Stomiidae","Gonostomatidae", "Sternoptychidae") ~"Stomiiformes", family %in%c("Platytroctidae","Alepocephalidae") ~"Alepocephaliformes", family =="Anoplogastridae"~"Trachichthyiformes", family %in%c("Paralepididae","Evermannellidae","Lestidiidae") ~"Aulopiformes", family =="Bathylagidae"~"Argentiniformes", family %in%c("Derichthyidae","Serrivomeridae") ~"Anguilliformes", family =="Eurypharyngidae"~"Saccopharyngiformes", family =="Zoarcidae"~"Perciformes", ) )
## Summary of the assemblages * species data.frame ----asb_sp_fish_summ <- mFD::asb.sp.summary(asb_sp_w = depth_fish_biomass)asb_sp_fish_occ <- asb_sp_fish_summ$"asb_sp_occ"htmltools::tagList(DT::datatable(asb_sp_fish_occ))
2.2 Computing distances between species based on functional traits
We have non-continuous traits so we use the Gower distance(metric = “gower”) as this method allows traits weighting.
2.3 Building functional spaces and chosing the best one
2.3.1 Computing several multimensional functional spaces and assessing their quality
mFD evaluates the quality of PCoA-based multidimensional spaces according to the deviation between trait-based distances and distances in the functional space (extension of Maire et al. (2015) framework).
2.3.3 Testing the correlation between functional axes and traits
Code
sp_faxes_coord_fish <- fspaces_quality_fish$"details_fspaces"$"sp_pc_coord"# As we have 26 traits we have to split the df to see correlation between functional axes and traits # first set ----fish_traits_1 <- fish_traits%>%select(1:9)fish_tr_faxes <- mFD::traits.faxes.cor(sp_tr = fish_traits_1, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes$"tr_faxes_stat"[which(fish_tr_faxes$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
1 eye_size PC1 Linear Model r2 0.161 0.0093
2 eye_size PC2 Linear Model r2 0.183 0.0052
3 eye_size PC3 Linear Model r2 0.133 0.0192
4 eye_size PC4 Linear Model r2 0.139 0.0163
5 orbital_length PC1 Linear Model r2 0.410 0.0000
10 gill_outflow PC2 Linear Model r2 0.511 0.0000
13 oral_gape_surface PC1 Linear Model r2 0.139 0.0165
14 oral_gape_surface PC2 Linear Model r2 0.513 0.0000
18 oral_gape_shape PC2 Linear Model r2 0.166 0.0083
24 oral_gape_position PC4 Linear Model r2 0.221 0.0019
26 lower_jaw_length PC2 Linear Model r2 0.699 0.0000
29 head_length PC1 Linear Model r2 0.325 0.0001
30 head_length PC2 Linear Model r2 0.376 0.0000
34 body_depth PC2 Linear Model r2 0.352 0.0000
35 body_depth PC3 Linear Model r2 0.203 0.0031
Code
## Plot ----fish_tr_faxes$"tr_faxes_plot"
Code
# second set ----fish_traits_2 <- fish_traits%>%select(10:18)fish_tr_faxes_2 <- mFD::traits.faxes.cor(sp_tr = fish_traits_2, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_2$"tr_faxes_stat"[which(fish_tr_faxes_2$"tr_faxes_stat"$"p.value"<0.05), ]
trait axis test stat value p.value
2 pectoral_fin_position PC2 Linear Model r2 0.175 0.0065
5 pectoral_fin_insertion PC1 Linear Model r2 0.279 0.0004
6 pectoral_fin_insertion PC2 Linear Model r2 0.388 0.0000
9 transversal_shape PC1 Linear Model r2 0.123 0.0246
11 transversal_shape PC3 Linear Model r2 0.220 0.0020
14 caudal_throttle_width PC2 Linear Model r2 0.410 0.0000
19 dorsal_fin_insertion PC3 Linear Model r2 0.184 0.0051
23 eye_position PC3 Linear Model r2 0.204 0.0030
32 ventral_photophores PC4 Kruskal-Wallis eta2 0.590 0.0000
33 gland_head PC1 Kruskal-Wallis eta2 0.245 0.0011
Code
## Plot ----fish_tr_faxes_2$"tr_faxes_plot"
Code
# third set ----fish_traits_3 <- fish_traits%>%select(19:26)fish_tr_faxes_3 <- mFD::traits.faxes.cor(sp_tr = fish_traits_3, sp_faxes_coord = sp_faxes_coord_fish[ , c("PC1", "PC2", "PC3", "PC4")], plot = T)## Print traits with significant effect ----fish_tr_faxes_3$"tr_faxes_stat"[which(fish_tr_faxes_3$"tr_faxes_stat"$"p.value"<0.05), ]
Standard Effect Size (SES): to eliminate the influence of species richness on the functional diversity indices (Mouchet et al., 2010). Measures the deviation from the random expectation in standard deviation units
null model frequency: Randomize community data matrix abundances (here biomasss) within species (maintains species occurrence frequency)
FRic Functional Richness: the proportion of functional space filled by species of the studied assemblage, i.e. the volume inside the convex-hull shaping species. To compute FRic the number of species must be at least higher than the number of functional axis + 1.
FDis Functional Dispersion: the biomass weighted deviation of species traits values from the center of the functional space filled by the assemblage i.e. the biomass-weighted mean distance to the biomass-weighted mean trait values of the assemblage.
FDiv Functional Divergence: the proportion of the biomass supported by the species with the most extreme functional traits i.e. the ones located close to the edge of the convex-hull filled by the assemblage.
FEve Functional Evenness: the regularity of biomass distribution in the functional space using the Minimum Spanning Tree linking all species present in the assemblage.
calcul des indices avec package FD
correction distance matrix (obtenue avec Gower) = “none” or “lingoes” transformation avec “sqrt” ne permet pas de transformer la matrice en “euclidean”
randomisation des traits des espèces par couche de profondeur
Code
# # Model parameters ----# # Number of simulations# n_simulations <- 9# # # correction method to use when the distance matrix cannot be represented in a Euclidean space# corr_method <- "lingoes"# "none"# # # Depth layers# depth_layers <- rownames(depth_fish_biomass)# # # List of indices to be calculated# indices <- c("FRic", "FDis", "FDiv", "FEve")# # # Matrices for storing observed and simulated results# dbFD_result_obs <- matrix(NA, nrow = length(depth_layers), ncol = length(indices),# dimnames = list(depth_layers, indices))# # dbFD_result_sim <- array(NA, dim = c(length(depth_layers), n_simulations, length(indices)),# dimnames = list(depth_layers, paste0("Sim.", 1:n_simulations), indices))# # # Function to randomize traits while preserving factor levels----# randomize_traits <- function(traits_mat) {# randomized_mat <- traits_mat# # Randomisation of columns# for (trait in colnames(traits_mat)) {# randomized_mat[, trait] <- sample(traits_mat[, trait], replace = FALSE)# }# return(randomized_mat)# }# # # Loop on each depth layer ----# for (layer in depth_layers) {# # # Select the species present in the selected depth layer# species_in_layer <- colnames(depth_fish_biomass)[depth_fish_biomass[layer, ] > 0]# # # Filter lines (traits) and remove constant columns# traits_layer <- fish_traits[species_in_layer, , drop = FALSE] %>%# select(where(~ length(unique(.)) > 1))# # # If no variable lines remain, we move on to the next layer# if (ncol(traits_layer) == 0) next# # # Correct formatting of biomass_layer with depth layer in rownames# # Converts biomass data to matrix, which is required by FD::dbFD()# biomass_layer <- matrix(as.numeric(depth_fish_biomass[layer, species_in_layer, drop = FALSE]),# nrow = 1, dimnames = list(layer, species_in_layer))# # # Calculation of observed indices with FD::dbFD# dbFD_result <- FD::dbFD(x = traits_layer, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,# corr = corr_method)# # # Storage of observed values# dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic# dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis# dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve# dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv# # # Simulations for each layer# for (sim in 1:n_simulations) {# randomized_traits <- randomize_traits(traits_layer)# # dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun = TRUE, m = 4, a = biomass_layer, messages=F,# corr = corr_method)# # # Storage of simulated values# dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic# dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis# dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve# dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv# }# }# # # Calculation of the SES (Standardized Effect Size) ----# SES_results <- data.frame(depth_layer = rownames(dbFD_result_obs))# # for (index in indices) {# meanNullFD <- rowMeans(dbFD_result_sim[, , index], na.rm = TRUE)# sdNullFD <- apply(dbFD_result_sim[, , index], 1, sd, na.rm = TRUE)# SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD# }# # # Plot ----# results_df <- tidyr::pivot_longer(# SES_results,# cols = -depth_layer,# names_to = "index",# values_to = "SES_FD"# )# # results_df$depth_layer <- factor(# results_df$depth_layer,# levels = c(# "Epipelagic",# "Upper mesopelagic",# "Lower mesopelagic",# "Bathypelagic"# )# )# # results_df$index <- factor(# results_df$index,# levels = c("FRic",# "FDis",# "FDiv",# "FEve"),# labels = c(# "Functional richness",# "Functional dispersion",# "Functional divergence",# "Functional evenness"# )# )# # ggplot(results_df, aes(x = depth_layer, y = SES_FD, color = depth_layer)) +# geom_point(size = 3) +# facet_wrap(~index) +# geom_hline(yintercept = c(1.96, -1.96), linetype = "dashed", color = "gray40", size=0.8) +# scale_color_manual(values = c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +# labs(x = "", y = "Standard Effect Size (SES)") +# theme_light() +# theme(axis.text.x = element_blank(),# strip.text.x = element_text(size = 14, color = "black"),# strip.background = element_rect(fill = "white"),# axis.title = element_text(size = 13),# axis.text = element_text(size = 13)) +# guides(col="none", fill="none")# # ggsave("SES_dbFD.png", path = "figures", dpi = 700, height = 5, width = 7)
3.2. Two levels (stations and depth layer)
Code
station_sp <-rbind(data_biomass_2002_2019, data_biomass_2021_2022) %>%as.data.frame() %>%left_join(metadata) %>%select(species, biomass_sp, volume_filtered, station) %>%# Divide biomass by the volume filtered at each trawl (g.m3)mutate(biomass_cpu = (biomass_sp / volume_filtered) *1000) %>%select(species, biomass_cpu, station) %>%replace(is.na(.), 0) %>%group_by(species, station) %>%mutate(biomass =sum(biomass_cpu)) %>%select(-biomass_cpu) %>%distinct() %>% tidyr::pivot_wider(names_from = species, values_from = biomass) %>%replace(is.na(.), 0) %>%arrange(station) %>%filter(!station %in%c("H0411", "L0731", "L0736")) %>% tibble::column_to_rownames(var ="station") %>%select(order(colnames(.))) %>%as.matrix()# Model parameters ----n_simulations <-9# Correction method for non-Euclidean distancescorr_method <-"lingoes"depth_layers <-rownames(station_sp) # Indices to calculateindices <-c("FRic", "FDis", "FDiv", "FEve")# Function to randomize traits while preserving factor levels----randomize_traits <-function(traits_mat) { randomized_mat <- traits_mat# Randomisation of columnsfor (trait incolnames(traits_mat)) { randomized_mat[, trait] <-sample(traits_mat[, trait], replace =FALSE) }return(randomized_mat)}# Data definition: list with depth and station matrices ----analysis_levels <-list(depth_layer = depth_fish_biomass,station = station_sp)# Store results ----all_results <-list()dbFD_result_sim_list <-list()# Function to calculate FD indices and SES ----calculate_FD_and_SES <-function(data_matrix, traits_data, n_simulations, indices, corr_method) { layers <-rownames(data_matrix) dbFD_result_obs <-matrix(NA, nrow =length(layers), ncol =length(indices), dimnames =list(layers, indices)) dbFD_result_sim <-array(NA, dim =c(length(layers), n_simulations, length(indices)),dimnames =list(layers, paste0("Sim.", 1:n_simulations), indices))# Loop through each layerfor (layer in layers) { species_in_layer <-colnames(data_matrix)[data_matrix[layer, ] >0] traits_layer <- traits_data[species_in_layer, , drop =FALSE] %>%select(where(~length(unique(.)) >1)) %>%droplevels()if (ncol(traits_layer) ==0) next biomass_layer <-matrix(as.numeric(data_matrix[layer, species_in_layer, drop =FALSE]),nrow =1, dimnames =list(layer, species_in_layer))# Calculate FD indices for observed data dbFD_result <- FD::dbFD(x = traits_layer, w.abun =TRUE, m =4, a = biomass_layer,messages =FALSE, corr = corr_method)# Verify structure of dbFD_result indices_names <-names(dbFD_result) # Extract valid indices valid_indices <-intersect(indices, indices_names) dbFD_result_obs[layer, valid_indices] <- dbFD_result[valid_indices]# Random simulationsfor (sim in1:n_simulations) { randomized_traits <-randomize_traits(traits_layer) dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun =TRUE, m =4, a = biomass_layer,messages =FALSE, corr = corr_method)# Verify structure of dbFD_sim_result indices_names_sim <-names(dbFD_sim_result) valid_indices_sim <-intersect(indices, indices_names_sim)# Store simulation results dbFD_result_sim[layer, sim, valid_indices_sim] <- dbFD_sim_result[valid_indices_sim] } }# Calculate SES (Standardized Effect Size) ---- SES_results <-data.frame(name =rownames(dbFD_result_obs))for (index in indices) { meanNullFD <-rowMeans(dbFD_result_sim[, , index], na.rm =TRUE) sdNullFD <-apply(dbFD_result_sim[, , index], 1, sd, na.rm =TRUE) SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD }return(list(dbFD_result_obs = dbFD_result_obs, dbFD_result_sim = dbFD_result_sim, SES_results = SES_results))}# Loop through analysis levelsfor (level innames(analysis_levels)) { data_matrix <- analysis_levels[[level]] layers <-rownames(data_matrix)# Matrices for observed and simulated results dbFD_result_obs <-matrix(NA, nrow =length(layers), ncol =length(indices),dimnames =list(layers, indices)) dbFD_result_sim <-array(NA, dim =c(length(layers), n_simulations, length(indices)),dimnames =list(layers, paste0("Sim.", 1:n_simulations), indices))# Loop through each layerfor (layer in layers) { species_in_layer <-colnames(data_matrix)[data_matrix[layer, ] >0] traits_layer <- fish_traits[species_in_layer, , drop =FALSE] %>%select(where(~length(unique(.)) >1)) %>%droplevels()if (ncol(traits_layer) ==0) next biomass_layer <-matrix(as.numeric(data_matrix[layer, species_in_layer, drop =FALSE]),nrow =1, dimnames =list(layer, species_in_layer))# Calculate FD indices for the observed data dbFD_result <- FD::dbFD(x = traits_layer, w.abun =TRUE, m =4, a = biomass_layer,messages =FALSE, corr = corr_method)# Store observed results dbFD_result_obs[layer, "FRic"] <- dbFD_result$FRic dbFD_result_obs[layer, "FDis"] <- dbFD_result$FDis dbFD_result_obs[layer, "FEve"] <- dbFD_result$FEve dbFD_result_obs[layer, "FDiv"] <- dbFD_result$FDiv# Perform random simulationsfor (sim in1:n_simulations) { randomized_traits <-randomize_traits(traits_layer) dbFD_sim_result <- FD::dbFD(x = randomized_traits, w.abun =TRUE, m =4, a = biomass_layer,messages =FALSE, corr = corr_method)# Store simulated results dbFD_result_sim[layer, sim, "FRic"] <- dbFD_sim_result$FRic dbFD_result_sim[layer, sim, "FDis"] <- dbFD_sim_result$FDis dbFD_result_sim[layer, sim, "FEve"] <- dbFD_sim_result$FEve dbFD_result_sim[layer, sim, "FDiv"] <- dbFD_sim_result$FDiv } }# Save simulated results for normality tests dbFD_result_sim_list[[level]] <- dbFD_result_sim# Calculate SES (Standardized Effect Size) SES_results <-data.frame(level = level, name =rownames(dbFD_result_obs))for (index in indices) { meanNullFD <-rowMeans(dbFD_result_sim[, , index], na.rm =TRUE) sdNullFD <-apply(dbFD_result_sim[, , index], 1, sd, na.rm =TRUE) SES_results[[index]] <- (dbFD_result_obs[, index] - meanNullFD) / sdNullFD } all_results[[level]] <- SES_results}# Normality and symmetry tests ----library(nortest)library(moments)# Normality and symmetry resultsnormality_results <-list()for (level innames(all_results)) { SES_data <- all_results[[level]] simulated_values <- dbFD_result_sim_list[[level]] normality_test <-data.frame(level = level, name = SES_data$name)for (index in indices) {# Normality test (Shapiro-Wilk) normality_test[[paste0(index, "_shapiro_p")]] <-apply(simulated_values[, , index], 1, function(x) shapiro.test(x)$p.value)# Skewness normality_test[[paste0(index, "_skewness")]] <-apply(simulated_values[, , index], 1, skewness) } normality_results[[level]] <- normality_test}# Stockage final des résultatsnormality_results_combined <-bind_rows(normality_results$station, normality_results$depth_layer, .id ="level")# Combine resultsSES_results_depth_layer <- all_results$depth_layer %>% tidyr::pivot_longer(cols =-c(name, level), names_to ="index", values_to ="SES_FD") %>%mutate(depth_layer=name)SES_results_station <- all_results$station %>%rename(station=name) %>% tidyr::pivot_longer(cols =-c(station, level), names_to ="index", values_to ="SES_FD") %>%inner_join(metadata %>%select(station, depth), by ="station") %>%mutate(depth_layer =case_when(between(depth, 0, 174) ~"Epipelagic",between(depth, 175, 699) ~"Upper mesopelagic",between(depth, 700, 999) ~"Lower mesopelagic",between(depth, 1000, 2000) ~"Bathypelagic" ) )# Combine datasetsSES_results_combined <-bind_rows(SES_results_station, SES_results_depth_layer)SES_results_combined$depth_layer <-factor( SES_results_combined$depth_layer,levels =c("Epipelagic","Upper mesopelagic","Lower mesopelagic","Bathypelagic" ))SES_results_combined$index <-factor( SES_results_combined$index,levels =c("FRic", "FDis", "FDiv", "FEve"),labels =c("Functional richness","Functional dispersion","Functional divergence","Functional evenness" ))# Plot ----ggplot(SES_results_combined, aes(x = depth_layer, y = SES_FD)) +geom_boxplot(data = SES_results_combined %>%filter(level =="station"),aes(color = depth_layer, fill = depth_layer),alpha =0.08, width =0.5, outlier.shape =NA) +geom_point(data = SES_results_combined %>%filter(level =="depth_layer"),aes(color = depth_layer),size =5, shape =18) +geom_jitter(data = SES_results_combined %>%filter(level =="station"),aes(color = depth_layer),size =1, width =0.2, alpha =0.35) +geom_hline(yintercept =c(1.96, -1.96), linetype ="dashed", color ="gray40", size =0.8) +scale_color_manual(values =c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +scale_fill_manual(values =c("#FEA520", "#D62246", "#6255B4", "#3C685A")) +facet_wrap(~ index, labeller =labeller(index =c(FRic ="Functional richness",FDis ="Functional dispersion",FDiv ="Functional divergence",FEve ="Functional evenness" ))) +labs(x =NULL, y ="Standard Effect Size (SES)") +theme_light() +theme(axis.text.x =element_blank(),strip.text.x =element_text(size =14, color ="black"),strip.background =element_rect(fill ="white"),axis.title =element_text(size =13),axis.text =element_text(size =13) ) +guides(col ="none", fill ="none")